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A method for performing wave packet simulations in dipole 
fields is presented. Starting from a Hamiltonian with non 
commuting terms, a gauge transformation leads to a new 
Hamiltonian which allows to calculate explicitly the evolu- 
tion operator. In this new gauge, the dipole field is fully 
included in the vector potential. The method of Goldberg, 
Schwartz and Schey based on the Caley form of the evolution 
operator is then generalized, and the resulting scheme is ap- 
plied to describe a switching device. The device is composed 
of a well region separated from a free region by a potential 
barrier, such that one bound state and one quasi level are 
present. For a particle initially in the ground state, transi- 
tions to the quasi level are induced by applying a dipole field, 
and the wave function can subsequently tunnel in the free re- 
gion. The probability to tunnel in the free region exhibits a 
plateaux structure as the wave function is emitted by "bursts" 
after each Rabi oscillation has been completed. 



The pioneering work of Goldberg, Schey and Schwartz 
[0 on numerical simulations of wave packet evolution 
through potential barriers and wells continues to have a 
deep impact on the teaching community. However, when 
this work first appeared, the technological means to ob- 
serve and characterize experimentally such phenomenon 
were absent. One dimensional quantum mechanics was 
then considered on a rather academic level. With the 
recent advances in the fabrication of semiconductor nano 
structures, one now has the means to design such struc- 
tures using molecular beam epitaxy techniques, etc... 
The motivation for research in this field now comes from 
the hope that quantum devices may one day be used 
as elementary building blocks for the next generation of 
computers. In the following, I will describe a method 
which allows to simulate the evolution of a wave packet 
in a potential landscape, in the presence of a dipole po- 
tential. The method will be applied to describe a switch- 
ing device consisting of a well, separated by a free region 
by a potential barrier. Starting with a particle which 
is trapped in the well, it is possible to induce a current 
in the free region by appling a micro wave field which 
activates the particle to a quasi level in the well. The 
tunneling current can be tuned by adjusting the excita- 
tion frequency, the amplitude of the external field and 
the characteristics of the barrier. 

While simulations of wave packet transmission through 
oscillating barriers have been achieved in the past pfl , lit- 



tle is known for the more realistic situation where the ap- 
plied potential depends linearly on the distance (electric 
field constant in space). This is relevant when a micro 
wave field is applied to a semiconductor heterostructure. 
Our starting point is the Hamiltonian H = Ho + 4>{x, t) 
with the free Hamiltonian: 



H 



2m dx 2 
and the dipole potential: 



V(x) 



(x, t) = ex cos(wi) 



(1) 



(2) 



with e the coupling strength and uj the driving frequency. 
For simplicity, I will consider static potentials V(x) which 
are constant in a set of different regions separated by po- 
tential steps. For more general situations, the potential 
V{x) can be decomposed as a series of infinitesimal steps 
[Q. Given an initial state ifi(x, 0), the wave function at 
later times is obtained by application of the evolution 
operator: 



ip(x,t) = U(t,0)ip(x,0) 

From standard Quantum Mechanics textbooks 
evolution operator is known to be: 



U(t,t ) = Texp 



dt'(H + <t>(x,t')) 



to 



(3) 
the 



(4) 



where the symbol T denotes the time ordering procedure. 
In past calculations, the assumption that the time vary- 
ing field <j)(t) be a constant in a given region of space 
allowed to drop the time ordering factor: Hq commutes 
with 4>(t) in this particular case. For the Hamiltonian of 
Eqs. (0) and (Q), the time ordering has to be kept because 
the kinetic energy does not commute with the position 
operator x. How can a numerical algorithm similar to 
that of Goldberg, Schey and Schwartz (GSS) Q be im- 
plemented? A "brute force" solution to this problem is to 
choose the discretization in time sufficiently small so that 
only the lowest order terms (0(H)) in the exponential of 
Eq. (f|) are important: the time ordering does not af- 
fect the approximate evolution operator for infinitesimal 
times: 



U(t + 5t,t) 



t+St 



dt'{H + cb{x,t)) (5) 



This ceases to be true when second order terms are taken 
into account. In addition to this blemish, another prob- 
lem arises: in situations where a large spatial region is 
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required, such as in the case where a (small) barrier or 
well region which specify the characteristic energies of 
the problem, is coupled to a (large) continuum region, 
the expectation value < x > of the position operator 
may take large values. This puts further restrictions on 
the magnitude of the time step for the simulation. 

Nevertheless, a more elegant method, which exploits 
the gauge invariance property of the Hamiltonian, does 
not suffer from the same constraints. The Hamiltonian 
of Eqs. (0) and (g) is written in a gauge where the vector 
potential A = and the scalar potential is (j>(x,t). An 
alternative choice of gauge can be obtained as: 



A' 



A- TpX 
ox 

x 1 d 

c at 



(6) 
(7) 



Now, I choose the "new" gauge so that the scalar poten- 
tial (jjl = 0. This then yields: 



X(x,t) 



■ sin(wt) 



A'(x,t) = — sin(urf) 



(8) 
(9) 



In this new gauge, the Hamiltonian is written as 

2 

I +V(x) 

k l UX C J 

n 2 d 2 



2m V i ox c 



+ i C ^ A' ® + V(a 
2m dx 2 mc dx 



t A 



12 



(10) 



2mc 2 

The second line follows from the fact that A' commutes 
with the momentum operator. In the new gauge, the 
wave function ip'(x,t) is related to the old one: 

ip'(x,t) = ij){x,t)e- iex{x > t)/hc (11) 

Note that the Hamiltonian H' does not suffer from the 
same setbacks as H: in each constant region of V(x) all 
terms in the Hamiltonian H' commute with each other, 
and consequently, the time ordering factor in the expres- 
sion for the evolution operator U'(t,to) can be dropped 
out: 



U'(t,to)=exp 



±J*dHH'{H) 



(12) 



exp 



t~t 
-i— — H 



g(Mo) d 

m dx 



i9(t,t ) 



with 



g(t,to) = — ^(cos(wt) — cos(wio)) 



0(Mo) 



2hmc 2 j t 



dt'A' 2 (t') 



(13) 
(14) 



Note that the last term in the exponential contributes 
only a time dependent phase factor to the evolution. 
Therefore, 8(t, to) will be omitted in the following. 



I now generalize the numerical scheme of GGS to the 
case where a vector potential is present. The elementary 
time step evolution for the wave function t/j' is taken in 
the so-called Caley form: 



-1 l_ Z-l TTl 

v V ; 1 + j-StH 1 v K ' 



(15) 



or, alternatively, in the "old" gauge: 



M + St) = e ie x(*+*>/ Rc - €^L e - iex ^' hc i)(t) (16) 

vy ' 1 + j-StH' 

which has the advantage of being exact to order (St) 2 
and unitary. Taking the convention h = 1, m = 1/2, and 
choosing the discretized variables x = jSx (j = 0, 1, ... J), 
t = nSt (n = 1,2...), (i>(x,t) = ^), (g(t) = 5tg n ), 
(V(x) = Vj), this is rewritten as: 

^+1(1 + i*!L Vj) + i\-\-^+t + 2^" +1 - ifpt) 

where A = 2Sx 2 /St. Introducing the quantity: 

f2« = - iSxg n ) + ^(2 + Sx 2 V 3 + iX) 

+^_ 1 (-l + iSxg n ) (18) 

Eq. (|l7| ) then takes the form: 



4>*+l(l + iSxg n ) + ^ +1 {-2 - Sx 2 V 3 + iX) 
+^+l(l-i6xg n ) = fl] 



(19) 



The above equation can be solved with the ansatz jjj : 



^"i 1 = e "^; +1 + fj 1 



(20) 



Which yields expressions for the quantities e™ and /": 



1 - iSxg n 1 2 + Sx 2 Vj — iX 
1 + iSxg n e"_ x 1 + iSxg n 

1 - iSxg n 



1 + iSxg n e^_-y 1 + iSxg rl 



-2lb 



a (21) 



(22) 



Boundary conditions for the above quantities are now 
needed. These are obtained from the boundary condition 
on the wave function (i/jq = ip™ = for all n): 



+1 _ 2 + 5x 2 V 1 -iX +1 _Jll 



1 + iSxg n 
which in turn implies 



1 + iSxg n 



(23) 
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fll 

Jl 



2 + fe 2 Ui - iX 
1 + iSxg n 

n't 



1 + iSxg„ 



(24) 
(25) 



Using Eqs. (pl]a— pT]b), one can now obtain e™ and /" 
for j = 2, .., J — 1. Similarly, the boundary condition for 



i> n ,t\ yields: 



(26) 



One can now update the wave function for j = J — 2, 
using the ansatz of Eq. (BO) : 



(27) 



which completes the numerical scheme. 

I now apply this method to a specific example. It has 
long been known || for a two level system that if a har- 
monic perturbation is applied such that the excitation 
frequency matches approximately the spacing between 
the two levels, the probability of occupation of a given 
level undergoes slow oscillations. Such oscillations are 
known in optics as Rabi oscillations. If the system is 
started in the ground state, and the dipole potential is 
switched on at t = 0, the system undergoes transitions to 
the excited state, but will ultimately return in the ground 
state after a time Tr = 2nh/\Fi2\ (resonant case), where 
F\2 is the matrix element of the harmonic perturbation. 
These slow oscillations will be exploited to construct a 
switch. 

Consider the potential landscape of Fig. pi a poten- 
tial well is separated by a "free" region by a potential 
barrier. The well bottom lies below the zero energy, and 
is adjusted so that there is only one bound state in the 
well. For energies E > 0, states extend from the free re- 
gion to the well region. If the barrier height was infinite 
a succession of bound states would be present in the well. 
However, due to the finite width of the barrier separating 
the well region from the free region, the excited states be- 
come quasi-levels, which have a finite lifetime in the well. 
In the numerical calculations which follow, we have ad- 
justed the height of the barrier such that only one quasi 
level is present among ~ 100 "continuum states" . By 
applying a dipole field on this system, with a frequency 
which corresponds to the spacing between the ground 
state and a quasi level, the wave function of a particle, 
initially in the ground state, will make transitions to the 
quasi level, and consequently leak out in the free region if 
the barrier is not too thick. The current generated in the 
free region then depends on the amplitude and frequency 
of the driving field, as well as the characteristics of the 
barrier. 

Typical values for the parameters of the potential of 
Fig. |l| are a = 4, b = 2, c = 200, V = 1, W = 3. The 
bound and excited states wave functions are specified by 
the connection formulas for the wave functions at each 



boundary, and the corresponding energies of these states 
are determined numerically. To determine the energy of 
the quasi levels, the integrated density in the well region 
for all states with E > is calculated, and the level for 
which this quantity is a maximum is selected. Alterna- 
tively, all matrix elements of the position operator be- 
tween the ground and excited states are computed, and 
the level for which the probability of transition is a max- 
imum is selected. In practice, these two procedures give 
the same result. Once the spacing between the quasi level 
energy Eq and the ground state energy Eq is known, the 
driving frequency u> is chosen to correspond to an exact 
resonance to — Eq — Eq, or alternatively, one can include 
a finite mismatch Slo. At t = 0, the particle is taken to 
be in the ground state, and at each time step, the inte- 
grated density pi{t) = J^t b +C dx p(x, t) in the free region 
is computed, as well as the overlap | < G\ip(t) > | 2 with 
the ground state. The time step has to be chosen small 
compared to the period of the external field: here, we 
choose St = 0.0125 x (2tt/u;). 

In Fig. ||a, pi(t) is plotted for field amplitude e = 0.1. 
The barrier width (b = 2) and height (W = 3) are cho- 
sen to be large enough that the "escape time" of the 
wave function is large compared to other characteristic 
times of the problem. Moreover, the driving frequency 
has been chosen to be close to the resonance condition 
(Slo/lu = 0.0001), which is smaller than the spacing be- 
tween "continuum" levels (E > 0), in order to check 
agreement with the two level approximation. The in- 
tegrated density exhibits steps or plateaux, which al- 
low to identify the Rabi frequency lor. Superposed to 
the plateaux structure are small amplitude oscillations 
which period corresponds to the driving frequency. As 
the simulation is started, transitions to the quasi lev- 
els and neighboring levels start occurring, but after a 
period Tr — 2w /lor, the contribution of the wave func- 
tion which remained trapped in the well has returned 
for the most part in the ground state. This is illus- 
trated in fig. ||b : indeed, aside from a slow decay as- 
sociated with the transparency of the barrier, the over- 
lap with the ground state | < G\ip(t) > | 2 oscillates 
with period which matches exactly the plateaux struc- 
ture. After the first plateau, pi(t) picks up again as 
the next oscillation returns the trapped wave function 
to the excited states. Upon doubling/halving the cou- 
pling strength, the period of the oscillations is twice as 
small/large, confirming the fact that the Rabi oscillation 
frequency scales linearly with the field amplitude if the 
resonance condition is met. The "measured" Rabi fre- 
quency ujft ~ 0.01 is in reasonable agreement with the 
two level result hujR = \ < G\x\l > |e ~ 0.04. To es- 
timate the magnitude of the matrix element < G|x|l > 
between the ground state and the quasi level, the param- 
eters of an infinite well were chosen. 

The evolution of the wave function for the same pa- 
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rameters is depicted in Figs. ||a-c for times t = 1500, 
t = 3000, and t = 4500, where the probability den- 
sity is plotted as a function of position. The well and 
barrier region lies on the left hand side of the pictures. 
The Rabi frequency corresponds roughly to a period 
~ 1500, and the leakage current is small, which ex- 
plains by most of the wave function remains in the well. 
At time t = 1500 (Fig. ||a), a portion of the wave func- 
tion has been transmitted in the free region, as a Rabi 
oscillation with the quasi level has been completed. At 
t = 3000 (Fig. ||b), the system has undergone two Rabi 
oscillations and another wave packet escapes from the 
well, and similarly for t = 4500 (Fig. ||c) after three 
oscillations. The wave function is therefore emitted by 
"bursts" ' out of the well every time a Rabi oscillation is 
completed. 

In summary, a method has been described, which simu- 
lates the quantum evolution of a particle in a square bar- 
rier/well potential, in the presence of a dipole field. Us- 
ing a gauge transformation, a Hamiltonian which terms 
commute with each other is obtained, which in turns al- 
lows to write an exact expression for the evolution opera- 
tor. Using a generalization of the finite difference scheme 
of (GSS), the time evolution is obtained. This method 
is particularly useful to model activation/tunneling pro- 
cesses for nano structures exposed to an external mi- 
crowave field. As an example, it is possible to use the 
two state slow (Rabi) oscillations of a two level system 
to build a switching device. The device is composed of a 
well region separated from a continuum region by a bar- 
rier. The activation of the ground state to a quasi level 
using microwave frequencies allows to generate a current 
in the free region in a controllable manner. More de- 
tails on the characteristics of this device will be provided 
elsewhere ||. 
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FIG. 1. Potential landscape considered: a well region is 
separated by a potential barrier from the "free region". The 
well depth and barrier height are adjusted so that there is only 
one bound state in the well, and there is only one quasi-level 
below the barrier 

FIG. 2. a) Integrated density in the "free" region as a func- 
tion of time, for e = 0.1. b) Overlap with the ground state as 
a function of time, with the same input parameters. 

FIG. 3. Probability density as a function of position. The 
well and the barrier are placed on the left hand side of the 
pictures, a) t = 1500 b) t = 3000 c) t = 4500 (t is measured 
in number of time steps St — 0.0125 x (2tt/ui). 
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